/*****************************************************************************************
 *                                                                                       *
 * OpenSpace                                                                             *
 *                                                                                       *
 * Copyright (c) 2014-2020                                                               *
 *                                                                                       *
 * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
 * software and associated documentation files (the "Software"), to deal in the Software *
 * without restriction, including without limitation the rights to use, copy, modify,    *
 * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
 * permit persons to whom the Software is furnished to do so, subject to the following   *
 * conditions:                                                                           *
 *                                                                                       *
 * The above copyright notice and this permission notice shall be included in all copies *
 * or substantial portions of the Software.                                              *
 *                                                                                       *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
 * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
 * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
 * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
 * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
 * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
 ****************************************************************************************/

/*****************************************************************************************
 * Modified parts of the code (4D texture mechanism) from Eric Bruneton is used in the   *
 * following code.                                                                       *
 ****************************************************************************************/

 /**
 * Precomputed Atmospheric Scattering
 * Copyright (c) 2008 INRIA
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the copyright holders nor the names of its
 *    contributors may be used to endorse or promote products derived from
 *    this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
 * THE POSSIBILITY OF SUCH DAMAGE.
 */

#version 460  core

; // preprocessor add semicolon to isolate error messages
#line 1 1 // C:\os\clean\OpenSpace-0.16.0\shaders\floatoperations.glsl
/*****************************************************************************************
 *                                                                                       *
 * OpenSpace                                                                             *
 *                                                                                       *
 * Copyright (c) 2014-2020                                                               *
 *                                                                                       *
 * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
 * software and associated documentation files (the "Software"), to deal in the Software *
 * without restriction, including without limitation the rights to use, copy, modify,    *
 * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
 * permit persons to whom the Software is furnished to do so, subject to the following   *
 * conditions:                                                                           *
 *                                                                                       *
 * The above copyright notice and this permission notice shall be included in all copies *
 * or substantial portions of the Software.                                              *
 *                                                                                       *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
 * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
 * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
 * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
 * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
 * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
 ****************************************************************************************/

#ifndef _FLOATOPERATIONS_GLSL_
#define _FLOATOPERATIONS_GLSL_

/**
 * Convert a positive floating point distance [0, 10^27]
 * (size of observable universe)
 * to a float in the range [-1, 1], suitable for depth buffer storage.
 * Note: This needs to be a monotonic function, so that the value can
 * still be used for depth comparison.
 */
float normalizeFloat(float inpt) {
    if (inpt > 1.0) {
        return inpt / pow(10, 30);
    } else {
        return inpt - 1.0;
    }
}

float denormalizeFloat(float inpt) {
    if (inpt < 0.0) {
        return inpt + 1.0;
    } else {
        return inpt * pow(10, 30);
    }
}

/**
 * Compute the length of a vector.
 * Supporting huge vectors, where the square of any of the components is too large to be
 * represented as a float.
 */
float safeLength(vec4 v) {
    float m = max(max(max(abs(v.x), abs(v.y)), abs(v.z)), abs(v.w));
    if (m > 0.f) {
        return length(v / m) * m;
    } else {
        return 0.f;
    }
}

float safeLength(vec3 v) {
    float m = max(max(abs(v.x), abs(v.y)), abs(v.z));
    if (m > 0.f) {
        return length(v / m) * m;
    } else {
        return 0.f;
    }
}

float safeLength(vec2 v) {
    float m = max(abs(v.x), abs(v.y));
    if (m > 0.f) {
        return length(v / m) * m;
    } else {
        return 0.f;
    }
}

/**
 * Normalize a vector
 * Supporting huge vectors, where the square of any of the components is too large to be
 * represent as a float.
 */
vec3 safeNormalize(vec3 v) {
    float m = max(max(abs(v.x), abs(v.y)), abs(v.z));
    return normalize(v / m);
}

#endif // _FLOATOPERATIONS_GLSL_
; // preprocessor add semicolon to isolate error messages
#line 63 0 // C:\os\clean\OpenSpace-0.16.0\modules\atmosphere\shaders\atmosphere_deferred_fs.glsl

; // preprocessor add semicolon to isolate error messages
#line 1 2 // C:\os\clean\OpenSpace-0.16.0\modules\atmosphere\shaders\atmosphere_common.glsl
/*****************************************************************************************
 *                                                                                       *
 * OpenSpace                                                                             *
 *                                                                                       *
 * Copyright (c) 2014-2020                                                               *
 *                                                                                       *
 * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
 * software and associated documentation files (the "Software"), to deal in the Software *
 * without restriction, including without limitation the rights to use, copy, modify,    *
 * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
 * permit persons to whom the Software is furnished to do so, subject to the following   *
 * conditions:                                                                           *
 *                                                                                       *
 * The above copyright notice and this permission notice shall be included in all copies *
 * or substantial portions of the Software.                                              *
 *                                                                                       *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
 * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
 * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
 * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
 * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
 * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
 ****************************************************************************************/

/*****************************************************************************************
 * Modified parts of the code (4D texture mechanism, analytical transmittance etc)       *
 * from Eric Bruneton is used in the following code.                                     *
 ****************************************************************************************/

 /**
 * Precomputed Atmospheric Scattering
 * Copyright (c) 2008 INRIA
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the copyright holders nor the names of its
 *    contributors may be used to endorse or promote products derived from
 *    this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
 * THE POSSIBILITY OF SUCH DAMAGE.
 */


// Atmosphere Rendering Parameters
uniform float Rg;
uniform float Rt;
uniform float AverageGroundReflectance;
uniform float groundRadianceEmittion;
uniform float HR;
uniform vec3 betaRayleigh;
uniform float HO;
uniform vec3 betaOzoneExtinction;
uniform float HM;
uniform vec3 betaMieScattering;
uniform vec3 betaMieExtinction;
uniform float mieG;
uniform float sunRadiance;

uniform bool ozoneLayerEnabled;

uniform int TRANSMITTANCE_W;
uniform int TRANSMITTANCE_H;
uniform int SKY_W;
uniform int SKY_H;
uniform int OTHER_TEXTURES_W;
uniform int OTHER_TEXTURES_H;
uniform int SAMPLES_R;
uniform int SAMPLES_MU;
uniform int SAMPLES_MU_S;
uniform int SAMPLES_NU;

const float ATM_EPSILON = 1.0;

// Integration steps
const int TRANSMITTANCE_STEPS = 500;
const int INSCATTER_INTEGRAL_SAMPLES = 50;
const int IRRADIANCE_INTEGRAL_SAMPLES = 32;
const int INSCATTER_SPHERICAL_INTEGRAL_SAMPLES = 16;

const float M_PI = 3.141592657;

uniform sampler2D transmittanceTexture;

float Rg2    = Rg * Rg;
float Rt2    = Rt * Rt;
float H      = sqrt(Rt2 - Rg2);
float H2     = Rt2 - Rg2;
float invSamplesMu = 1.0f / float(SAMPLES_MU);
float invSamplesR = 1.0f / float(SAMPLES_R);
float invSamplesMuS = 1.0f / float(SAMPLES_MU_S);
float invSamplesNu = 1.0f / float(SAMPLES_NU);
float RtMinusRg = float(Rt - Rg);
float invRtMinusRg = 1.0f / RtMinusRg;

float opticalDepth(const float localH, const float r, const float mu, const float d) {
  float invH = 1.0/localH;
  float a    = sqrt((0.5 * invH)*r);
  vec2 a01   = a*vec2(mu, mu + d / r);
  vec2 a01s  = sign(a01);
  vec2 a01sq = a01*a01;
  float x    = a01s.y > a01s.x ? exp(a01sq.x) : 0.0;
  vec2 y     = a01s / (2.3193*abs(a01) + sqrt(1.52*a01sq + 4.0)) * vec2(1.0, exp(-d*invH*(d/(2.0*r)+mu)));
  return sqrt((6.2831*H)*r) * exp((Rg-r)*invH) * (x + dot(y, vec2(1.0, -1.0)));
}

vec3 analyticTransmittance(const float r, const float mu, const float d) {
  if (ozoneLayerEnabled) {
    return exp(-betaRayleigh * opticalDepth(HR, r, mu, d) -
             betaOzoneExtinction * (0.0000006) * opticalDepth(HO, r, mu, d) -
             betaMieExtinction * opticalDepth(HM, r, mu, d));
  }
  else {
    return exp(-betaRayleigh * opticalDepth(HR, r, mu, d) -
               betaMieExtinction * opticalDepth(HM, r, mu, d));
  }
}

vec3 irradiance(sampler2D sampler, const float r, const float muSun) {
  float u_r     = (r - Rg) * invRtMinusRg;
  float u_muSun = (muSun + 0.2) / (1.0 + 0.2);
  return texture(sampler, vec2(u_muSun, u_r)).rgb;
}


//================================================//
//=============== General Functions ==============//
//================================================//
// In the following shaders r (altitude) is the length of vector/position x in the
// atmosphere (or on the top of it when considering an observer in space),
// where the light is comming from the opposite direction of the view direction,
// here the vector v or viewDirection.
// Rg is the planet radius and Rt the atmosphere radius.

//--- Calculate the distance of the ray starting at x (height r)
// until the planet's ground or top of atmosphere. ---
// r := || vec(x) || e [0, Rt]
// mu := cosine of the zeith angle of vec(v). Or mu = (vec(x) * vec(v))/r
float rayDistance(const float r, const float mu) {
  // The light ray starting at the observer in/on the atmosphere can
  // have to possible end points: the top of the atmosphere or the
  // planet ground. So the shortest path is the one we are looking for,
  // otherwise we may be passing through the ground.
  
  // cosine law
  float atmRadiusEps = Rt + ATM_EPSILON;
  float rayDistanceAtmosphere = -r * mu +
    sqrt(r * r * (mu * mu - 1.0f) + atmRadiusEps * atmRadiusEps);
  float delta = r * r * (mu * mu - 1.0f) + Rg * Rg;

  // Ray may be hitting ground
  if (delta >= 0.0f) {
    float rayDistanceGround = -r * mu - sqrt(delta);
    if (rayDistanceGround >= 0.0f) {
      return min(rayDistanceAtmosphere, rayDistanceGround);
    }
  }
  return rayDistanceAtmosphere;
}

//-- Given the window's fragment coordinates, for a defined
// viewport, gives back the interpolated r e [Rg, Rt] and
// mu e [-1, 1] --
// r := height of starting point vect(x)
// mu := cosine of the zeith angle of vec(v). Or mu = (vec(x) * vec(v))/r
void unmappingRAndMu(out float r, out float mu) {
  float u_mu  = gl_FragCoord.x / float(TRANSMITTANCE_W);
  float u_r   = gl_FragCoord.y / float(TRANSMITTANCE_H);
  
  // In the paper u_r^2 = (r^2-Rg^2)/(Rt^2-Rg^2)
  // So, extracting r from u_r in the above equation:
  //r  = sqrt( Rg * Rg + (u_r * u_r) * (Rt * Rt - Rg * Rg) );
  r  = Rg + (u_r * u_r) * RtMinusRg;
  
  // In the paper the Bruneton suggest mu = dot(v,x)/||x|| with ||v|| = 1.0
  // Later he proposes u_mu = (1-exp(-3mu-0.6))/(1-exp(-3.6))
  // But the below one is better. See Colliene.
  // One must remember that mu is defined from 0 to PI/2 + epsillon.
  mu = -0.15f + tan(1.5f * u_mu) / tan(1.5f) * (1.0f + 0.15f);
}

//-- Given the windows's fragment coordinates, for a defined view port,
// gives back the interpolated r e [Rg, Rt] and muSun e [-1, 1] --
// r := height of starting point vect(x)
// muSun := cosine of the zeith angle of vec(s). Or muSun = (vec(s) * vec(v))
void unmappingRAndMuSun(out float r, out float muSun) {
  // See Bruneton and Colliene to understand the mapping.
  muSun = -0.2f + (gl_FragCoord.x - 0.5f) / (float(OTHER_TEXTURES_W) - 1.0f) * (1.0f + 0.2f);
  //r  = Rg + (gl_FragCoord.y - 0.5f) / (float(OTHER_TEXTURES_H) - 1.0f) * (Rt - Rg);
  r  = Rg + (gl_FragCoord.y - 0.5f) / (float(OTHER_TEXTURES_H) ) * RtMinusRg;
}

//-- Given the windows's fragment coordinates, for a defined view port,
// gives back the interpolated r e [Rg, Rt] and muSun e [-1, 1] for the
// Irradiance deltaE texture table --
// r := height of starting point vect(x)
// muSun := cosine of the zeith angle of vec(s). Or muSun = (vec(s) * vec(v))
void unmappingRAndMuSunIrradiance(out float r, out float muSun) {
  // See Bruneton and Colliene to understand the mapping.
  muSun = -0.2f + (gl_FragCoord.x - 0.5f) / (float(SKY_W) - 1.0f) * (1.0f + 0.2f);
  r  = Rg + (gl_FragCoord.y - 0.5f) / (float(SKY_H) - 1.0f) * RtMinusRg;
}

//-- Given the windows's fragment coordinates, for a defined view port,
// gives back the interpolated r e [Rg, Rt] and mu, muSun amd nu e [-1, 1] --
// r := height of starting point vect(x)
// mu := cosine of the zeith angle of vec(v). Or mu = (vec(x) * vec(v))/r
// muSun := cosine of the zeith angle of vec(s). Or muSun = (vec(s) * vec(v))
// nu := cosone of the angle between vec(s) and vec(v)
// dhdH := it is a vec4. dhdH.x stores the dminT := Rt - r, dhdH.y stores the dH value (see paper),
// dhdH.z stores dminG := r - Rg and dhdH.w stores dh (see paper).
void unmappingMuMuSunNu(const float r, vec4 dhdH, out float mu, out float muSun, out float nu) {
  // Window coordinates of pixel (uncentering also)
  float fragmentX = gl_FragCoord.x - 0.5f;
  float fragmentY = gl_FragCoord.y - 0.5f;

  // Pre-calculations
  //float Rg2 = Rg * Rg;
  //float Rt2 = Rt * Rt;
  float r2  = r * r;
  
  float halfSAMPLE_MU = float(SAMPLES_MU) / 2.0f;
  // If the (vec(x) dot vec(v))/r is negative, i.e.,
  // the light ray has great probability to touch
  // the ground, we obtain mu considering the geometry
  // of the ground
  if (fragmentY < halfSAMPLE_MU) {
    float ud = 1.0f - (fragmentY / (halfSAMPLE_MU - 1.0f));
    float d  = min(max(dhdH.z, ud * dhdH.w), dhdH.w * 0.999);
    // cosine law: Rg^2 = r^2 + d^2 - 2rdcos(pi-theta)
    // where cosine(theta) = mu
    mu = (Rg2 - r2 - d * d) / (2.0 * r * d);
    // We can't handle a ray inside the planet, i.e.,
    // when r ~ Rg, so we check against it.
    // If that is the case, we approximate to
    // a ray touching the ground.
    // cosine(pi-theta) = dh/r = sqrt(r^2-Rg^2)
    // cosine(theta) = - sqrt(1 - Rg^2/r^2)
    mu = min(mu, -sqrt(1.0 - (Rg2 / r2)) - 0.001);
  }
  // The light ray is touching the atmosphere and
  // not the ground
  else {
    float d = (fragmentY - halfSAMPLE_MU) / (halfSAMPLE_MU - 1.0f);
    d = min(max(dhdH.x, d * dhdH.y), dhdH.y * 0.999);
    // cosine law: Rt^2 = r^2 + d^2 - 2rdcos(pi-theta)
    // whre cosine(theta) = mu
    mu = (Rt2 - r2 - d * d) / (2.0f * r * d);
  }
  
  float modValueMuSun = mod(fragmentX, float(SAMPLES_MU_S)) / (float(SAMPLES_MU_S) - 1.0f);
  // The following mapping is different from the paper. See Colliene for an details.
  muSun = tan((2.0f * modValueMuSun - 1.0f + 0.26f) * 1.1f) / tan(1.26f * 1.1f);
  nu = -1.0f + floor(fragmentX / float(SAMPLES_MU_S)) / (float(SAMPLES_NU) - 1.0f) * 2.0f;
}


//-- Function to access the transmittance texture. Given r
// and mu, returns the transmittance of a ray starting at vec(x),
// height r, and direction vec(v), mu, and length until it hits
// the ground or the top of atmosphere. --
// r := height of starting point vect(x)
// mu := cosine of the zeith angle of vec(v). Or mu = (vec(x) * vec(v))/r
vec3 transmittanceLUT(const float r, const float mu) {
  // Given the position x (here the altitude r) and the view
  // angle v (here the cosine(v)= mu), we map this
  float u_r  = sqrt((r - Rg) * invRtMinusRg);
  //float u_r  = sqrt((r*r - Rg*Rg) / (Rt*Rt - Rg*Rg));
  // See Colliene to understand the different mapping.
  float u_mu = atan((mu + 0.15f) / (1.0f + 0.15f) * tan(1.5f)) / 1.5f;
  
  return texture(transmittanceTexture, vec2(u_mu, u_r)).rgb;
}

// -- Given a position r and direction mu, calculates de transmittance
// along the ray with length d. This function uses the propriety
// of Transmittance: T(a,b) = TableT(a,v)/TableT(b, v) --
// r := height of starting point vect(x)
// mu := cosine of the zeith angle of vec(v). Or mu = (vec(x) * vec(v))/r
vec3 transmittance(const float r, const float mu, const float d) {
  // Here we use the transmittance property: T(x,v) = T(x,d)*T(d,v)
  // to, given a distance d, calculates that transmittance along
  // that distance starting in x (hight r): T(x,d) = T(x,v)/T(d,v).
  //
  // From cosine law: c^2 = a^2 + b^2 - 2*a*b*cos(ab)
  float ri = sqrt(d * d  + r * r + 2.0 * r * d * mu);
  // mu_i = (vec(d) dot vec(v)) / r_i
  //      = ((vec(x) + vec(d-x)) dot vec(v))/ r_i
  //      = (r*mu + d) / r_i
  float mui = (d + r * mu) / ri;

  // It's important to remember that we calculate the Transmittance
  // table only for zenith angles between 0 and pi/2+episilon.
  // Then, if mu < 0.0, we just need to invert the view direction
  // and the start and end points between them, i.e., if
  // x --> x0, then x0-->x.
  // Also, let's use the property: T(a,c) = T(a,b)*T(b,c)
  // Because T(a,c) and T(b,c) are already in the table T,
  // T(a,b) = T(a,c)/T(b,c).
  if (mu > 0.0f) {
    return min(transmittanceLUT(r, mu) /
               transmittanceLUT(ri, mui), 1.0f);
  }
  else {
    return min(transmittanceLUT(ri, -mui) /
               transmittanceLUT(r, -mu), 1.0f);
  }
}

// -- Calculates Rayleigh phase function given the
// scattering cosine angle mu --
// mu := cosine of the zeith angle of vec(v). Or mu = (vec(x) * vec(v))/r
float rayleighPhaseFunction(const float mu) {
    //return (3.0f / (16.0f * M_PI)) * (1.0f + mu * mu);
    return 0.0596831036 * (1.0f + mu * mu);
}

// -- Calculates Mie phase function given the
// scattering cosine angle mu --
// mu := cosine of the zeith angle of vec(v). Or mu = (vec(x) * vec(v))/r
float miePhaseFunction(const float mu) {
  //return (3.0f / (8.0f * M_PI)) *
  //      ( ( (1.0f - (mieG * mieG) ) * (1.0f + mu * mu) ) /
  //      ( (2.0f + mieG * mieG) *
  //        pow(1.0f + mieG * mieG - 2.0f * mieG * mu, 3.0f/2.0f) ) );
  // return 1.5f * 1.0f / (4.0f * M_PI) * (1.0f - mieG * mieG) *
  //   pow(1.0f + (mieG * mieG) - 2.0f * mieG * mu, -3.0f/2.0f) * (1.0f + mu * mu) / (2.0f + mieG*mieG);

  float mieG2 = mieG * mieG;
  return 0.1193662072 * (1.0f - mieG2) *
    pow(1.0f + mieG2 - 2.0f * mieG * mu, -1.5f) * (1.0f + mu * mu) / (2.0f + mieG2);
}

// -- Given the height rm view-zenith angle (cosine) mu,
// sun-zenith angle (cosine) muSun and the angle (cosine)
// between the vec(s) and vec(v), nu, we access the 3D textures
// and interpolate between them (r) to find the value for the
// 4D texture. --
// r := height of starting point vect(x)
// mu := cosine of the zeith angle of vec(v). Or mu = (vec(x) * vec(v))/r
// muSun := cosine of the zeith angle of vec(s). Or muSun = (vec(s) * vec(v))
// nu := cosine of the angle between vec(s) and vec(v)
vec4 texture4D(sampler3D table, const float r, const float mu,
                const float muSun, const float nu)
{
  //float Rg2    = Rg * Rg;
  //float Rt2    = Rt * Rt;
  float r2     = r * r;
  //float H      = sqrt(Rt2 - Rg2);
  float rho    = sqrt(r2 - Rg2);
  float rmu    = r * mu;
  float delta  = rmu * rmu - r2 + Rg2;
  //float invSamplesMu = 1.0f / float(SAMPLES_MU);
  //float invSamplesR = 1.0f / float(SAMPLES_R);
  //float invSamplesMuS = 1.0f / float(SAMPLES_MU_S);
  //float invSamplesNu = 1.0f / float(SAMPLES_NU);
  // vec4 cst     = rmu < 0.0f && delta > 0.0f ?
  //   vec4(1.0f, 0.0f, 0.0f, 0.5f - 0.5f / float(SAMPLES_MU)) :
  //   vec4(-1.0f, H * H, H, 0.5f + 0.5f / float(SAMPLES_MU));

  vec4 cst     = rmu < 0.0f && delta > 0.0f ?
    vec4(1.0f, 0.0f, 0.0f, 0.5f - 0.5f * invSamplesMu) :
    vec4(-1.0f, H2, H, 0.5f + 0.5f * invSamplesMu);

  //float u_r    = 0.5f / float(SAMPLES_R) + rho / H * (1.0f - 1.0f / float(SAMPLES_R));
  float u_r    = 0.5f * invSamplesR + rho / H * (1.0f - invSamplesR);
  //float u_mu   = cst.w + (rmu * cst.x + sqrt(delta + cst.y)) / (rho + cst.z) * (0.5f - 1.0f / float(SAMPLES_MU));
  float u_mu   = cst.w + (rmu * cst.x + sqrt(delta + cst.y)) / (rho + cst.z) * (0.5f - invSamplesMu);
  // float u_mu_s = 0.5f / float(SAMPLES_MU_S) +
  //   (atan(max(muSun, -0.1975) * tan(1.26f * 1.1f)) / 1.1f + (1.0f - 0.26f)) * 0.5f * (1.0f - 1.0f / float(SAMPLES_MU_S));
  float u_mu_s = 0.5f * invSamplesMuS +
    (atan(max(muSun, -0.1975) * tan(1.386f)) * 0.9090909090909090 + (0.74f)) * 0.5f * (1.0f - invSamplesMuS);
  float lerp = (nu + 1.0f) / 2.0f * (float(SAMPLES_NU) - 1.0f);
  float u_nu = floor(lerp);
  lerp = lerp - u_nu;

  // return texture(table, vec3((u_nu + u_mu_s) / float(SAMPLES_NU), u_mu, u_r)) * (1.0f - lerp) +
  //   texture(table, vec3((u_nu + u_mu_s + 1.0f) / float(SAMPLES_NU), u_mu, u_r)) * lerp;
  
  return texture(table, vec3((u_nu + u_mu_s) * invSamplesNu, u_mu, u_r)) * (1.0f - lerp) +
    texture(table, vec3((u_nu + u_mu_s + 1.0f) * invSamplesNu, u_mu, u_r)) * lerp;
}

// -- Given the irradiance texture table, the cosine of zenith sun vector
// and the height of the observer (ray's stating point x), calculates the
// mapping for u_r and u_muSun and returns the value in the LUT. --
// lut   := OpenGL texture2D sampler (the irradiance texture deltaE)
// muSun := cosine of the zeith angle of vec(s). Or muSun = (vec(s) * vec(v))
// r     := height of starting point vect(x)
vec3 irradianceLUT(sampler2D lut, const float muSun, const float r) {
  // See Bruneton paper and Coliene to understand the mapping
  float u_muSun = (muSun + 0.2f) / (1.0f + 0.2f);
  float u_r     = (r - Rg) * invRtMinusRg;
  return texture(lut, vec2(u_muSun, u_r)).rgb;
}
; // preprocessor add semicolon to isolate error messages
#line 65 0 // C:\os\clean\OpenSpace-0.16.0\modules\atmosphere\shaders\atmosphere_deferred_fs.glsl

out vec4 renderTarget;
in vec3 interpolatedNDCPos;
in vec2 texCoord;

uniform int nAaSamples;
uniform int cullAtmosphere;

uniform sampler2D irradianceTexture;
uniform sampler3D inscatterTexture;
uniform sampler2D mainPositionTexture;
uniform sampler2D mainNormalTexture;
uniform sampler2D mainColorTexture;

uniform dmat4 dInverseModelTransformMatrix;
uniform dmat4 dModelTransformMatrix;
uniform dmat4 dSGCTViewToWorldMatrix;
uniform dmat4 dSgctProjectionToModelTransformMatrix;

uniform dvec4 dCamPosObj;
uniform dvec3 sunDirectionObj;

/*******************************************************************************
 ***** ALL CALCULATIONS FOR ECLIPSE ARE IN METERS AND IN WORLD SPACE SYSTEM ****
 *******************************************************************************/
// JCC: Remove and use dictionary to
// decides the number of shadows
const uint numberOfShadows = 1;

struct ShadowRenderingStruct {
        double xu, xp;
        double rs, rc;
        dvec3 sourceCasterVec;
        dvec3 casterPositionVec;
        bool isShadowing;
};

// Eclipse shadow data
// JCC: Remove and use dictionary to
// decides the number of shadows
uniform ShadowRenderingStruct shadowDataArray[numberOfShadows];
uniform int shadows;
uniform bool hardShadows;

vec4 butterworthFunc(const float d, const float r, const float n) {
    return vec4(vec3(sqrt(r/(r + pow(d, 2*n)))), 1.0);
}

vec4 calcShadow(const ShadowRenderingStruct shadowInfoArray[numberOfShadows], const dvec3 position,
                const bool ground) {
    if (shadowInfoArray[0].isShadowing) {
        dvec3 pc = shadowInfoArray[0].casterPositionVec - position;
        dvec3 sc_norm = shadowInfoArray[0].sourceCasterVec;
        dvec3 pc_proj = dot(pc, sc_norm) * sc_norm;
        dvec3 d = pc - pc_proj;
        
        float length_d = float(length(d));
        double length_pc_proj = length(pc_proj);
        
        float r_p_pi = float(shadowInfoArray[0].rc * (length_pc_proj + shadowInfoArray[0].xp) / shadowInfoArray[0].xp);
        float r_u_pi = float(shadowInfoArray[0].rc * (shadowInfoArray[0].xu - length_pc_proj) / shadowInfoArray[0].xu);
        
        if ( length_d < r_u_pi ) { // umbra
            if (ground) {
              if (hardShadows) {
                  return vec4(0.2, 0.2, 0.2, 1.0);
              }
              else {
                  return butterworthFunc(length_d, r_u_pi, 2.0);
              }
            }
            else {
              if (hardShadows) {
                  return vec4(0.5, 0.5, 0.5, 1.0);
              }
              else {
                  //return vec4(vec3(length_d/r_p_pi), 1.0);
                  return butterworthFunc(length_d, r_u_pi, 2.0);
              }
            }
        }
        else if ( length_d < r_p_pi ) {// penumbra
            if (hardShadows) {
                return vec4(0.5, 0.5, 0.5, 1.0);
            }
            else {
                return vec4(vec3(length_d/r_p_pi), 1.0);
            }
        }
    }
     
    return vec4(1.0);
}

/*******************************************************************************
 ******* ALL CALCULATIONS FOR ATMOSPHERE ARE KM AND IN WORLD SPACE SYSTEM ******
 *******************************************************************************/

struct dRay {
    dvec4 origin;
    dvec4 direction;
};

/* Function to calculate the initial intersection of the eye (camera) ray
 * with the atmosphere.
 * In (all parameters in the same coordinate system and same units):
 * - planet position
 * - ray direction (normalized)
 * - eye position
 * - atmosphere radius
 * Out: true if an intersection happens, false otherwise
 * - inside: true if the ray origin is inside atmosphere, false otherwise
 * - offset: the initial intersection distance from eye position when
 *           the eye is outside the atmosphere
 * - maxLength : the second intersection distance from eye position when the
 *               eye is outside the atmosphere or the initial (and only)
 *               intersection of the ray with atmosphere when the eye position
 *               is inside atmosphere.
 */
bool dAtmosphereIntersection(const dvec3 planetPosition, const dRay ray, const double atmRadius,
                             out bool inside, out double offset, out double maxLength ) {
    dvec3  l  = planetPosition - ray.origin.xyz;
    double s  = dot(l, ray.direction.xyz);
    double l2 = dot(l, l);
    double r2 = atmRadius * atmRadius; // avoiding surface acne

    // Ray origin (eye position) is behind sphere
    if ((s < 0.0) && (l2 > r2)) {
        inside    = false;
        offset    = 0.0;
        maxLength = 0.0;
        return false;
    }

    double m2 = l2 - s*s;

    // Ray misses atmospere
    if (m2 > r2) {
        inside    = false;
        offset    = 0.0;
        maxLength = 0.0;
        return false;
    }

    // We already now the ray hits the atmosphere

    // If q = 0.0f, there is only one intersection
    double q = sqrt(r2 - m2);

    // If l2 < r2, the ray origin is inside the sphere
    if (l2 > r2) {
        inside    = false;
        offset    = s - q;
        maxLength = s + q;
    }
    else {
        inside    = true;
        offset    = 0.0;
        maxLength = s + q;
    }
    
    return true;
}

/*
 * Calculates Intersection Ray by walking through
 * all the graphic pipile transformations in the
 * opposite direction.
 * Instead of passing through all the pipeline,
 * it starts at NDC from the interpolated
 * positions from the screen quad.
 * This method avoids matrices multiplications
 * wherever is possible.
 */
void dCalculateRayRenderableGlobe(out dRay ray,
                                  out dvec4 planetPositionObjectCoords,
                                  out dvec4 cameraPositionInObject) {
    dvec4 clipCoords = dvec4(interpolatedNDCPos.xy, 1.0, 1.0);

    // Clip to Object Coords
    dvec4 objectCoords = dSgctProjectionToModelTransformMatrix * clipCoords;

    objectCoords /= objectCoords.w;
    
    // Planet Position in Object Space
    // JCC: Applying the inverse of the model transformation on the object postion in World
    // space results in imprecision.
    planetPositionObjectCoords = dvec4(0.0, 0.0, 0.0, 1.0);

    // Camera Position in Object Space (in meters)
    cameraPositionInObject = dCamPosObj;
    
    // ============================
    // ====== Building Ray ========
    // Ray in object space (in KM)
    ray.origin    = cameraPositionInObject * dvec4(0.001, 0.001, 0.001, 1.0);
    //ray.direction = dvec4(normalize(objectCoords.xyz - cameraPositionInObject.xyz), 0.0);
    ray.direction = dvec4(normalize((objectCoords.xyz * dvec3(0.001))- ray.origin.xyz), 0.0);
}

/*
 * Calculates the light scattering in the view direction comming from other
 * light rays scattered in the atmosphere.
 * Following the paper:  S[L]|x - T(x,xs) * S[L]|xs
 * The view direction here is the ray: x + tv, s is the sun direction,
 * r and mu the position and zenith cosine angle as in the paper.
 * Arguments:
 * x := camera position
 * t := ray displacement variable after calculating the intersection with the
 * atmosphere. It is the distance from the camera to the last intersection with
 * the atmosphere. If the ray hits the ground, t is updated to the correct value
 * v := view direction (ray's direction) (normalized)
 * s := Sun direction (normalized)
 * r := out of ||x|| inside atmosphere (or top of atmosphere)
 * mu := out of cosine of the zenith view angle
 * attenuation := out of transmittance T(x,x0). This will be used later when
 * calculating the reflectance R[L].
 */
vec3 inscatterRadiance(inout vec3 x, inout float t, inout float irradianceFactor,
                       const vec3 v, const vec3 s, out float r, out float mu,
                       out vec3 attenuation, const vec3 fragPosObj, out bool groundHit,
                       const double maxLength, const double pixelDepth,
                       const vec4 spaceColor, const float sunIntensity) {

    const float INTERPOLATION_EPS = 0.004f; // precision const from Brunetton

    vec3 radiance;
    
    r  = length(x);
    mu = dot(x, v) / r;

    float mu2           = mu * mu;
    float r2            = r * r;
    float nu            = dot(v, s);
    float muSun         = dot(x, s) / r;
    float rayleighPhase = rayleighPhaseFunction(nu);
    float miePhase      = miePhaseFunction(nu);
      
    // S[L](x,s,v)
    // I.e. the next line has the scattering light for the "infinite" ray passing
    // through the atmosphere. If this ray hits something inside the atmosphere,
    // we will subtract the attenuated scattering light from that path in the
    // current path.
    vec4 inscatterRadiance = max(texture4D(inscatterTexture, r, mu, muSun, nu), 0.0);
      
    // After removing the initial path from camera pos to top of atmosphere (for an
    // observer in the space) we test if the light ray is hitting the atmosphere
    vec3  x0     = fragPosObj;
    float r0     = length(fragPosObj);
    float invr0  = 1.0/r0;
    float muSun0 = dot(fragPosObj, s) * invr0;
    //vec3  x0     = x + float(pixelDepth) * v;
    float mu0    = dot(x0, v) * invr0;

    if ((pixelDepth > INTERPOLATION_EPS) && (pixelDepth < maxLength)) {
        t = float(pixelDepth);
        groundHit = true;
        
        // Transmittance from point r, direction mu, distance t
        // By Analytical calculation
        //attenuation = analyticTransmittance(r, mu, t);
        // JCC: change from analytical to LUT transmittance to avoid
        // acme on planet surface when looking from far away. (11/02/2017)
        attenuation = transmittance(r, mu, t);
        
        // Here we use the idea of S[L](a->b) = S[L](b->a), and get the S[L](x0, v, s)
        // Then we calculate S[L] = S[L]|x - T(x, x0)*S[L]|x0
        // The "infinite" ray hist something inside the atmosphere, so we need to remove
        // the unsused contribution to the final radiance.
        vec4 inscatterFromSurface = texture4D(inscatterTexture, r0, mu0, muSun0, nu);
        inscatterRadiance = max(inscatterRadiance - attenuation.rgbr * inscatterFromSurface, 0.0);

        // We set the irradianceFactor to 1.0 so the reflected irradiance will be considered
        // when calculating the reflected light on the ground.
        irradianceFactor = 1.0;
    }
    else {
        attenuation = analyticTransmittance(r, mu, t);
        //attenuation = transmittance(r, mu, t);
        groundHit = false;
    }

    // cos(PI-thetaH) = dist/r
    // cos(thetaH) = - dist/r
    // muHorizon = -sqrt(r^2-Rg^2)/r = -sqrt(1-(Rg/r)^2)
    float muHorizon = -sqrt(1.0f - (Rg2 / r2));

    // In order to avoid imprecision problems near horizon,
    // we interpolate between two points: above and below horizon
    if (abs(mu - muHorizon) < INTERPOLATION_EPS) {
        // We want an interpolation value close to 1/2, so the
        // contribution of each radiance value is almost the same
        // or it has a heavy weight if from above or below horizon
        float interpolationValue = ((mu - muHorizon) + INTERPOLATION_EPS) / (2.0f * INTERPOLATION_EPS);

        //float t2 = t * t;
              
        // Above Horizon
        mu  = muHorizon - INTERPOLATION_EPS;
        //r0  = sqrt(r * r + t * t + 2.0f * r * t * mu);
        // From cosine law where t = distance between x and x0
        // r0^2 = r^2 + t^2 - 2 * r * t * cos(PI-theta)
        //r0  = sqrt(r2 + t2 + 2.0f * r * t * mu);
        float halfCossineLaw1 = r2 + (t * t);
        float halfCossineLaw2 = 2.0f * r * t;
        r0  = sqrt(halfCossineLaw1 + halfCossineLaw2 * mu);
        
        float invr0 = 1.0/r0;
        // From the dot product: cos(theta0) = (x0 dot v)/(||ro||*||v||)
        // mu0 = ((x + t) dot v) / r0
        // mu0 = (x dot v + t dot v) / r0
        // mu0 = (r*mu + t) / r0
        mu0 = (r * mu + t) * invr0;
        
        vec4 inScatterAboveX  = texture4D(inscatterTexture, r, mu, muSun, nu);
        vec4 inScatterAboveXs = texture4D(inscatterTexture, r0, mu0, muSun0, nu);
        // Attention for the attenuation.r value applied to the S_Mie
        vec4 inScatterAbove = max(inScatterAboveX - attenuation.rgbr * inScatterAboveXs, 0.0f);

        // Below Horizon
        mu  = muHorizon + INTERPOLATION_EPS;
        //r0  = sqrt(r2 + t2 + 2.0f * r * t * mu);
        r0  = sqrt(halfCossineLaw1 + halfCossineLaw2 * mu);
        invr0 = 1.0/r0;
        
        mu0 = (r * mu + t) * invr0;
        
        vec4 inScatterBelowX  = texture4D(inscatterTexture, r, mu, muSun, nu);
        vec4 inScatterBelowXs = texture4D(inscatterTexture, r0, mu0, muSun0, nu);
        // Attention for the attenuation.r value applied to the S_Mie
        vec4 inScatterBelow = max(inScatterBelowX - attenuation.rgbr * inScatterBelowXs, 0.0);

        // Interpolate between above and below inScattering radiance
        inscatterRadiance = mix(inScatterAbove, inScatterBelow, interpolationValue);
    }

    // The w component of inscatterRadiance has stored the Cm,r value (Cm = Sm[L0])
    // So, we must reintroduce the Mie inscatter by the proximity rule as described in the
    // paper by Bruneton and Neyret in "Angular precision" paragraph:
      
    // Hermite interpolation between two values
    // This step is done because imprecision problems happen when the Sun is slightly below
    // the horizon. When this happens, we avoid the Mie scattering contribution.
    inscatterRadiance.w *= smoothstep(0.0f, 0.02f, muSun);
    vec3 inscatterMie    = inscatterRadiance.rgb * inscatterRadiance.a / max(inscatterRadiance.r, 1e-4) *
      (betaRayleigh.r / betaRayleigh);
      
    radiance = max(inscatterRadiance.rgb * rayleighPhase + inscatterMie * miePhase, 0.0f);
    
    // Finally we add the Lsun (all calculations are done with no Lsun so
    // we can change it on the fly with no precomputations)
    // return radiance * sunRadiance;
    vec3 finalScatteringRadiance = radiance * sunIntensity;
    
    if (groundHit) {
        return finalScatteringRadiance;
    }
    else {
        //return ((r-Rg) * invRtMinusRg)*spaceColor.rgb + finalScatteringRadiance;
        return spaceColor.rgb + finalScatteringRadiance;
        // return attenuation * spaceColor.rgb +
        // (vec3(1.0) - attenuation) * finalScatteringRadiance;
    }
}

/*
 * Calculates the light reflected in the view direction comming from other
 * light rays integrated over the hemispehre plus the direct light (L0) from Sun.
 * Following the paper: R[L]= R[L0]+R[L*]
 * The the ray is x + tv, v the view direction, s is the sun direction,
 * r and mu the position and zenith cosine angle as in the paper.
 * As for all calculations in the atmosphere, the center of the coordinate system
 * is the planet's center of coordiante system, i.e., the planet's position is (0,0,0).
 * Arguments:
 * x := camera position
 * t := ray displacement variable. Here, differently from the inScatter light calculation,
 * the position of the camera is already offset (on top of atmosphere) or inside
 * the atmosphere.
 * v := view direction (ray's direction) (normalized)
 * s := Sun direction (normalized)
 * r := ||x|| inside atmosphere (or top of atmosphere). r <= Rt here.
 * mu := cosine of the zenith view angle
 * attenuationXtoX0 := transmittance T(x,x0)
 */
vec3 groundColor(const vec3 x, const float t, const vec3 v, const vec3 s, const float r,
                 const float mu, const vec3 attenuationXtoX0, const vec4 groundColor,
                 const vec3 normal, const float irradianceFactor,
                 const float waterReflectance, const float sunIntensity)
{
    vec3 reflectedRadiance = vec3(0.0f);

    // First we obtain the ray's end point on the surface
    vec3  x0                 = x + t * v;
    float r0                 = length(x0);
    // Normal of intersection point.
    // Normal must be normalized.
    vec3  n                  = normal;
    //vec4 groundReflectance = groundColor * vec4(.37);
    vec4 groundReflectance   = groundColor *
        vec4(groundRadianceEmittion, groundRadianceEmittion, groundRadianceEmittion, 1.0f);

    // L0 is not included in the irradiance texture.
    // We first calculate the light attenuation from the top of the atmosphere
    // to x0.
    float dotNS = dot(n, s);
    float muSun = max(dotNS, 0.0f);

    // Is direct Sun light arriving at x0? If not, there is no direct light from Sun (shadowed)
    vec3  transmittanceL0     = muSun < -sqrt(1.0f - (Rg2 / (r0 * r0))) ?
                                vec3(0.0f) : transmittanceLUT(r0, muSun);
    // E[L*] at x0
    vec3  irradianceReflected = irradiance(irradianceTexture, r0, muSun) * irradianceFactor;

    // R[L0] + R[L*]
    // vec3 groundRadiance = (dotNS < -0.2f ? groundReflectance.rgb * 15 : groundReflectance.rgb) *
    //   (muSun * transmittanceL0 + irradianceReflected) * sunIntensity / M_PI;

    vec3 groundRadiance;
    vec3 RLStar = (muSun * transmittanceL0 + irradianceReflected) * sunIntensity / M_PI;
    if (dotNS < 0.05f) {
        groundRadiance = groundReflectance.rgb * mix(30.0f, 1.0f, smoothstep(-1.0f, 0.05f, dotNS)) * RLStar;
    }
    else {
        groundRadiance = groundReflectance.rgb * RLStar;
    }

    // Specular reflection from sun on oceans and rivers
    if ((waterReflectance > 0.1f) && /*(dotNS > -0.2f)*/(muSun > 0.0f)) {
        vec3  h         = normalize(s - v);
        // Fresnell Schlick's approximation
        float fresnel   = 0.02f + 0.98f * pow(1.0f - dot(-v, h), 5.0f);
        // Walter BRDF approximation
        float waterBrdf = fresnel * pow(max(dot(h, n), 0.0f), 150.0f);
        // Adding Fresnell and Water BRDFs approximation to the final surface color
        // (After adding the sunRadiance and the attenuation of the Sun through atmosphere)
        groundRadiance += waterReflectance * max(waterBrdf, 0.0f) * transmittanceL0 * sunIntensity;
    }
    //return groundRadiance;
    // Finally, we attenuate the surface Radiance from the the point x0 to the camera location.
    reflectedRadiance = attenuationXtoX0 * groundRadiance;
      
    // Returns reflectedRadiance = 0.0 if the ray doesn't hit the ground.
    return reflectedRadiance;
}

/*
 * Calculates the Sun color.
 * The the ray is x + tv, v the view direction, s is the sun direction,
 * r and mu the position and zenith cosine angle as in the paper.
 * As for all calculations in the atmosphere, the center of the coordinate system
 * is the planet's center of coordiante system, i.e., the planet's position is (0,0,0).
 * Arguments:
 * x := camera position
 * t := ray displacement variable. Here, differently from the inScatter light calculation,
 * the position of the camera is already offset (on top of atmosphere) or inside
 * the atmosphere.
 * v := view direction (ray's direction) (normalized)
 * s := Sun direction (normalized)
 * r := ||x|| inside atmosphere (or top of atmosphere). r <= Rt here.
 * mu := cosine of the zenith view angle
 * attenuation := transmittance T(x,x0)
 */
vec3 sunColor(const vec3 x, const float t, const vec3 v, const vec3 s, const float r,
              const float mu, const float irradianceFactor) {
    vec3 transmittance  = (r <= Rt) ? ( mu < -sqrt(1.0f - Rg2/(r*r)) ?
                          vec3(0.0f) : transmittanceLUT(r, mu)) : vec3(1.0f);
    // JCC: Change this function to a impostor texture with gaussian decay color weighted
    // by tge sunRadiance, transmittance and irradianceColor (11/03/2017)
    float sunFinalColor = smoothstep(cos(M_PI / 500.0f), cos(M_PI / 900.0f), dot(v, s)) *
                          sunRadiance * (1.0f - irradianceFactor);

    return transmittance * sunFinalColor;
}

void main() {
    if (cullAtmosphere == 0) {
        vec4 atmosphereFinalColor = vec4(0.0f);
        int nSamples = 1;
        
        // First we determine if the pixel is complex (different fragments on it)
        bool complex = false;
        vec4 oldColor, currentColor;
        vec4 colorTexture;
        
        colorTexture = texture(mainColorTexture, texCoord);
        
        // Color from G-Buffer
        vec4 color = colorTexture;
        // Ray in object space
        dRay ray;
        dvec4 planetPositionObjectCoords = dvec4(0.0);
        dvec4 cameraPositionInObject     = dvec4(0.0);
        
        // Get the ray from camera to atm in object space
        dCalculateRayRenderableGlobe(ray, planetPositionObjectCoords,
                                    cameraPositionInObject);
        
        bool  insideATM    = false;
        double offset      = 0.0;   // in Km
        double maxLength   = 0.0;   // in Km

        bool  intersectATM = false;

        intersectATM = dAtmosphereIntersection(planetPositionObjectCoords.xyz, ray,
                                                Rt - (ATM_EPSILON * 0.001), insideATM, offset, maxLength);
            
        if ( intersectATM ) {
            // Now we check is if the atmosphere is occluded, i.e., if the distance to the pixel
            // in the G-Buffer positions is less than the distance to the atmosphere then the atmosphere
            // is occluded
            // Fragments positions into G-Buffer are written in SGCT Eye Space (View plus Camera Rig Coords)
            // when using their positions later, one must convert them to the planet's coords
            
            // =======================
            // Get data from G-Buffer
            // =======================

            // Normal is stored in SGCT View Space and transformed to the current object space
            vec4 normalViewSpaceAndWaterReflectance = texture(mainNormalTexture, texCoord);
            dvec4 normalViewSpace = vec4(normalViewSpaceAndWaterReflectance.xyz, 0.0);
            dvec4 normalWorldSpace = dSGCTViewToWorldMatrix * normalViewSpace;
            vec4 normal = vec4(dInverseModelTransformMatrix * normalWorldSpace);
            normal.xyz = normalize(normal.xyz);
            normal.w = normalViewSpaceAndWaterReflectance.w;

            // Data in the mainPositionTexture are written in view space (view plus camera rig)
            vec4 position = texture(mainPositionTexture, texCoord);

            // OS Eye to World coords
            dvec4 positionWorldCoords = dSGCTViewToWorldMatrix * position;

            // World to Object (Normal and Position in meters)
            dvec4 positionObjectsCoords = dInverseModelTransformMatrix * positionWorldCoords;

            
            // Distance of the pixel in the gBuffer to the observer
            // JCC (12/12/2017): AMD distance function is buggy.
            //double pixelDepth = distance(cameraPositionInObject.xyz, positionObjectsCoords.xyz);
            double pixelDepth = length(cameraPositionInObject.xyz - positionObjectsCoords.xyz);
            
            // JCC (12/13/2017): Trick to remove floating error in texture.
            // We see a squared noise on planet's surface when seeing the planet
            // from far away.
            float dC = float(length(cameraPositionInObject.xyz));
            float x1 = 1e8;
            if (dC > x1) {
                pixelDepth     += 1000.0;
                float alpha     = 1000.0;
                float beta      = 1000000.0;
                float x2        = 1e9;
                float diffGreek = beta - alpha;
                float diffDist  = x2 - x1;
                float varA      = diffGreek/diffDist;
                float varB      = (alpha - varA * x1);
                pixelDepth     += double(varA * dC + varB);
            }

            // All calculations are done in Km:
            pixelDepth                *= 0.001;
            positionObjectsCoords.xyz *= 0.001;

            if (pixelDepth < offset) {
                // ATM Occluded - Something in fron of ATM.
                atmosphereFinalColor += color;
            }
            else {
                // Following paper nomenclature
                double t = offset;
                vec3 attenuation;

                // Moving observer from camera location to top atmosphere
                // If the observer is already inside the atm, offset = 0.0
                // and no changes at all.
                vec3  x  = vec3(ray.origin.xyz + t*ray.direction.xyz);
                float r  = 0.0f;//length(x);
                vec3  v  = vec3(ray.direction.xyz);
                float mu = 0.0f;//dot(x, v) / r;
                vec3  s  = vec3(sunDirectionObj);
                float tF = float(maxLength - t);

                // Because we may move the camera origin to the top of atmosphere
                // we also need to adjust the pixelDepth for tdCalculateRayRenderableGlobe' offset so the
                // next comparison with the planet's ground make sense:
                pixelDepth -= offset;
                
                dvec4 onATMPos           = dModelTransformMatrix * dvec4(x * 1000.0, 1.0);
                vec4 eclipseShadowATM    = calcShadow(shadowDataArray, onATMPos.xyz, false);
                float sunIntensityInscatter = sunRadiance * eclipseShadowATM.x;

                float irradianceFactor = 0.0;

                bool groundHit = false;
                vec3 inscatterColor = inscatterRadiance(x, tF, irradianceFactor, v,
                                                        s, r, mu, attenuation,
                                                        vec3(positionObjectsCoords.xyz),
                                                        groundHit, maxLength, pixelDepth,
                                                        color, sunIntensityInscatter);
                vec3 groundColorV = vec3(0.0);
                vec3 sunColorV = vec3(0.0);
                if (groundHit) {
                    vec4 eclipseShadowPlanet = calcShadow(shadowDataArray, positionWorldCoords.xyz, true);
                    float sunIntensityGround = sunRadiance * eclipseShadowPlanet.x;
                    groundColorV = groundColor(x, tF, v, s, r, mu, attenuation,
                                                color, normal.xyz, irradianceFactor,
                                                normal.a, sunIntensityGround);
                }
                else {
                    // In order to get better performance, we are not tracing
                    // multiple rays per pixel when the ray doesn't intersect
                    // the ground.
                    sunColorV = sunColor(x, tF, v, s, r, mu, irradianceFactor);
                }
                
                // Final Color of ATM plus terrain:
                vec4 finalRadiance = vec4(inscatterColor + groundColorV + sunColorV, 1.0);

                atmosphereFinalColor += finalRadiance;
            }
        }
        else { // no intersection
            // Buffer color
            atmosphereFinalColor += color;
        }
        

        renderTarget = atmosphereFinalColor;
    }
    else { // culling
        vec4 bColor = texture(mainColorTexture, texCoord);
        renderTarget = bColor;
    }
}

